Expansion of ID polarized superfiuids: The FFLO state reveals itself 



(N 

o 

(N 
P-h 



(Z3 1 
I 

^— > ■ 
+-» ■ 

C3 1 



o 
o 



>: 

in, 

ON' 

o: 



X 



Hong Lu 1 , L. 0. Baksmaty 1 , C. J. Bolech 2,1 , and Han Pu 1 
1 Department of Physics and Astronomy, Rice University, 6100 Main street, MS-61, Houston, TX 77005. 

2 Department of Physics, University of Cincinnati, 
345 Clifton Court, ML-11, Cincinnati, OH 45221-0011. 
(Dated: February 20, 2012) 

We study the expansion dynamics of a one dimensional polarized Fermi gas after sudden release 
from confinement using both the mean-field Bogoliubov-de Gennes and the numerically exact Time- 
Evolving Block Decimation methods. Our results show that experimentally observable spin density 
modulations directly related to the presence of a Fulde-Ferrel-Larkin-Ovchinnikov (FFLO) state 
develop during the expansion of the cloud, providing incontrovertible evidence of this long-sought 
state. 
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Since the introduction of the Bardeen-Cooper- 
Schrieffer (BCS) theory, physicists have speculated on 
the fate of the superconducting pairing correlation in the 
presence of a polarizing effect. This could arise from 
a mass imbalance of the pairing fermions such as in 
color superconductivity or in the vicinity of magnetic 
impurities within conventional superconductors. The 
FFLO (Fulde-Ferrel-Larkin-Ovchinnikov) 0-S1 proposal 
suggests that in such circumstances the condensation will 
occur from pairs with finite center-of-mass momenta. De- 
spite decades of work [U-El], this state has not been unam- 
biguously observed. Although recent experiments [7( in 
one dimension (ID) confirmed important aspects of the 
phase diagram [H, |f|, conclusive evidence of the FFLO 
phase was not obtained. We show here that during a non- 
equilibrium expansion, the polarized ID superfluid devel- 
ops strong signatures in the density profiles of the paring 
species which are a direct consequence of the FFLO crys- 
talline order and constitute incontrovertible evidence. 

We focus on a polarized degenerate Fermi gas con- 
fined to a ID harmonic trap. In general, according to 
the Mermin-Ho- Wagner theorem, a ID superfluid sys- 
tem cannot support superfluidity and would possess, at 
best, algebraically decaying long range order at zero tem- 
perature (T = 0). However for the finite systems that 
we study here, there is copious theoretical evidence that 
FFLO correlations occur and are fairly robust [lol- 17|. 



We also note that the experiments use not a single ID 
trap but a loosely coupled array which allows tuning of 
the inter-tube coupling and thus makes it possible to 
study the 3D to ID crossover physics jjj. Although a par- 
tially polarzied phase was observed through direct imag- 
ing in the experiment, it is quite clear from recent work 
that the FFLO correlations do not leave a detectable 
signature on the ground state density profiles. Thus the 
character of the partially polarized phase remains un- 
known. 

We consider a gas of N fermionic particles each of 
mass m with two spin projections labeled by a = (j, J,) 
confined to a cigar-shaped harmonic trap. Consistent 
with experimental reality 0, [l8l - l2^ |. we assume that 
the inter-particle interaction arises from a broad fesh- 



bach resonance and is amenable to exquisite control. In 
these systems, the ratio of the radial u> r and axial to z 
trapping frequencies which define the anisotropy of the 
trap A = U) r /oj z can be made so large that the Fermi 
energy Ep associated with the axial dynamics of the 
trap Nhuj z and the temperature kgT, are both much 
smaller than the energy level spacing of the radial con- 
finement fku r i.e., Nhuj z ,kBT << hu r 0. Due to ex- 
tremely rarerified nature of the gas, there are virtually 
no spin relaxation processes and the particles interact via 
s-wave scattering gir>S(z) . Furthermore, in addition to 
the total number TV, the total polarization of the cloud 
P = (JV t - Ni)/(N-f + iVj.) can also be varied through 
independent control of the number of particles in each 
spin projection N a . Formally this system is described by 
a Hamiltonian H = j dz (Hq + Hi) with : 



H (z) = 



h 2 d 2 

2m dz 2 



+ V~ trap (z) - p 



(1) 



where ip a (z) and [i a represent the fermionic field oper- 
ators and the chemical potential of atomic species with 
spin a and V^^z) — ^uj 2 z 2 . We define the Fermi en- 
ergy, radius, momentum and temperature as Ep = N, 
zp = y/2Ep, kp — \/2Ep and Tp = Ep and measure 
the relative strength of the interaction with the ratio (7) 
of the interaction (ej) and the kinetic (efc) energy densi- 
ties . In the limit of weak interaction ei ~ gmp(z) and 
e k - p 2 (z) yielding: 



7 = 9vd/p 



(2) 



Our calculations are done using two methods with dis- 
tinct but complementary advantages. First is the Time 
Evolving Block Decimation (TEBD) ^ (See Supple- 
mental Material at for details of methods) , an exact ap- 
proach that retains all important correlations. Second 
is the mean-field Bogoliubov-de Gennes (BdG) method, 
an effective theory approach which retains only the two 
point correlations and describes the spin densities p a {z) 
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Figure 1: (color online)The expansion of sample with TV = 40, P = 0.05 and <?id = —8.0. From left to right, each column 
represents snapshots of the expansion dynamics at t=Q.Q, 1.0, 1.7 (l/ui z ). Row 1 displays the density profiles. In each plot, we 
show p-f, pi and S = p-f — pi obtained from BdG calculation. Row 2 is the same as Row 1 except that the results are obtained 
from TEBD calculation. Row 3 shows the spin densities S(z) from the TEBD. Finally in Row 4 we plot the amplitude of the 
superfluid gap |A| from the BdG calculation. 



and the superfluid gap A(z) through quasi-particle wave- 
functions. The BdG has the advantage that, when cor- 
rect, it provides a clear picture of the dynamics of the 
pairing field A(z) = giD{tp^{z)ipi{z)) in direct associa- 
tion with the particle densities p a {z) = (z) t/v (z)}- 
However, although the BdG has been observed to give a 
very good description of ID samples at weak interaction 
[Io| . we do not expect this trend to extend from moder- 
ate to strong interaction. On the other hand the TEBD 
method provides a stringent check for the observed phe- 
nomena in the BdG approach. In both cases we work at 
T = and employ a canonical approach which fixes TV 
and P. 

To observe the FFLO state, experiments must ver- 
ify crystalline order in A(z) or alternatively, that the 
average center-of-mass momentum of the pairs (rife) is 
proportional to the separation of the Fermi surfaces 
(rife) oc fef — fcj,. In ID this relationship can be recast 
in terms of the spin density S(z) = p j [{z) — p\.{z) as 
(rife) oc 7r JV S(z)dz/ L, where L is the size of the par- 
tially polarized region. Recently a number of authors 
[111 , [la |24| have suggested the measurement of the pair 
momentum distribution function rife as the most promis- 
ing avenue to detecting the finite center-of-mass momen- 
tum q of the pairs. These suggestions are extrapolations 
from equilibrium studies where rik shows peaks at k = ±q 
in contrast to the peak at k = expected for regular 
BCS pairing. However, we are not aware of analyses of 
rife accounting for the interacting nature of the expansion 
dynamics and in particular how well this signal will be 
preserved. This is particularly important for ID given 



that 7 increases during expansion [see Eq. (|2|)]. In this 
study we explore the possibility of finding a signal di- 
rectly in real space. Our calculations reveal that: 

• Upon axial expansion, strong accents develop in the 
spin density profiles. 

• The position of these accents exactly coincide with 
the nodes in the pair correlation function and rep- 
resent prima facie evidence of FFLO correlations. 

• The strength of this signal increases with 7 and 
decreases with polarization, being strongest when 
the spin excitations are gapped. 

• The accents in the spin density move much more 
slowly than the edge of the cloud. 

In Fig. [1] dramatic accents in the spin densities are ob- 
served as the expansion of the cloud proceeds. Through 
a comparison of the density plots with the correspond- 
ing gap parameter |A(z)| (bottom row in Fig.[T]) one can 
make a key observation: The position and growth of the 
spin density accents respectively coincide with the nodes 
and amplification of\A(z)\. Furthermore, these spin den- 
sity accents (or the order parameter nodes) move much 
slower during the expansion compared to the edge of the 
whole cloud. 

To understand this phenomenon, it is helpful to first 
layout some broad features of the ground state utilizing 
the phase diagram for a homogeneous system together 
with the local density approximation (LDA) [H, 0, fl4j . 
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Figure 2: Density profiles, obtained from TEBD calculation, 
during the expansion of a sample with N = 40, P = 0.15 and 
gin = -8.0. 
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Figure 3: Pair momentum distribution at two different polar- 
ization for giD = —8 and N = 40. In each panel, we display 
rife for different times. Counting from the left, the curves cor- 
respond to t = 0, 0.47, 0.94 and 1.41, from top to bottom. 
In both cases the momentum peaks representing the FFLO 
state disappear from the plot during the expansion. 



Under LDA, the trapped system can be regarded as lo- 
cally homogeneous with chemical potential denned by: 
There are two regimes to be consid- 



in the pair momentum distribution rik defined by: 



A*( 2 ) = ,_ 
ered 0,113,111 



" Vtrap(z)- 



14j depending on whether P is smaller 
or larger than a critical polarization P c . For P < P c , 
we obtain an FFLO state at the center of the trap sur- 
rounded by fully paired BCS wings at the edges. Here the 
BdG calculation tells us that there is exactly one excess 
spin bound to each of the nodes of the order parameter 
and the FFLO state is analogous to a band insulator of 
the relative motion between the unpaired and paired par- 
ticles. The ground state represented in Fig. [T] is within 
this regime and density accents represent the localization 
of unpaired spins at the nodes of A. During the time 
of flight, the excess spins are kept pinned to the nodes 
of the order parameter and become more tightly bound. 
The dramatic effects observed occur when this localiza- 
tion couples with the average enhancement of | A | implied 
by an increasing 7 as the density drops during expansion 
[see Eq. fl5}]; a uniquely ID phenomenon. Henceforth we 
refer to these accents as node signatures. 

For P > P c , the FFLO state still remains at the cen- 
ter in the ground state, but the wings exclusively contain 
the majority spin component. In this regime, there are 
more excess spins than nodes of A, and consequently they 
are less tightly bound. Here we expect the node signa- 
tures to be less dramatic which is confirmed in Fig. [2] In 
particular, the spin accents near the edges are not well 
resolved. We can therefore conclude that the best place 
to observe the node signature is at P < P c where the sig- 
nal is enhanced by both a large separation of the nodes 
and greater contrast with the background density. We 
note that the value of P c increases with |<7id| implying a 
sizable observation window at strong interactions where 
experiments are conducted. 

At equilibrium the FFLO correlation appear as peaks 



- f f dzdz' e lk{z - z ' ] 0{z,z r ) , (3) 



where 0(z,z') = (ipt(z)tp^(z)ip^(z')ipf(z')) is the two- 
point correlation function. In Fig.[3]we observe the effects 
of interaction on this signature during the expansion. At 
sufficiently long time, nu no longer possesses peaks at 
finite momentum. 

One may wonder whether the node signatures can be 
observed in in situ density profiles of a trapped cloud 
with sufficiently large interaction strength. To answer 
this, we show in Fig. [4] the density profiles of a trapped 
system for giD = —8, —20 and —36. (Note that for the 
experiment reported in Ref. 0, <?id ~ — 50 for the cen- 
tral tube.) One can see that the modulation depth of 
the spin density of a trapped cloud is not very sensitive 
to giD- This is in sharp contrast to the BdG calculation 
where the spin density modulation is indeed enhanced 
as 7 is increased — an indication of the invalidity of the 
mean-field theory for strong interaction. In the exact cal- 
culation, the localization of excess spin at large |<7id| is 
counter-balanced by increased quantum fluctuations ne- 
glected in the mean-field theory. Therefore, the dramatic 
emergence of node signatures is a unique feature of ex- 
pansion dynamics. 

Finally, we address the question of the effect of the in- 
teraction strength in Fig. [5] where the spin densities in 
an expanding cloud are shown for two sets of interac- 
tion strength. Though the results from the strong and 
weak interaction are qualitatively similar, the spin ac- 
cents start to develop earlier for the case of smaller gn). 
This could play an important role in practice when finite 
lifetime of the system must be taken into account. 

In conclusion, we have investigated the expansion dy- 
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Figure 4: Ground state density profiles in trap, with AT = 40, 
P = 0.05 and for different interaction strengths gm. In each 
plot, we show pf, pi and S obtained from TEBD calculation. 
For clarity, the spin density S is magnified by a factor of 2. 
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Figure 5: Expansion profiles for two different samples with 
N = 40, P — 0.05 but at different interaction strengths gm- 
In each plot, we plot the TEBD result for S. In both cases, 
the modulation depth of the spin density first reduces and 
then strengthens during expansion. 



namics of polarized Fermi superfluid in ID using both 
the BdG and TEBD methods. Our results predict that 
strong spin density modulations which can be readily ob- 
served in experiment, emerge during expansion and pro- 
vide direct evidence of the FFLO state. Apart from the 
pair momentum distribution function described above, 
other methods 25] have been proposed in the literature 
to detect FFLO. However they all rely on interferomet- 
ric techniques requiring two fermionic superfhhds, one 
of them being the FFLO state. Our proposal, in con- 
trast, only requires the FFLO cloud itself and hence is 
significantly simpler. In a more general context, our work 
shows that quantum dynamics of low-dimensional atomic 
gases is highly non-trivial and deserve a more thorough 
study in the future. 
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METHODS - SUPPLEMENTARY MATERIAL 



This system of N = + N± is described by a Hamil- 
tonian H = J dz (Hq + Hj) with non-interacting (Hq) 
and interaction (Hi) energy densities given by: 



n 2 d 2 



H (z) = 

a 

Hi{z) = g 1D ip(z)ipl(z)il; l (z)ip t (z) 



(4) 



where 4 , a (z) represent the Fermionic field operators, m 
the mass and \i a the chemical potential of atomic species 
with spin a . The ID effective coupling constant giu < 
is expressed through a relationship with the 3D scattering 

Here a; is the 



length a 3D by [f]: g 1D 



2h a 3D 



mai(l-Aa 3Dd /ai) ' 

oscillator length and A = ^(f/2)/y / 2- We work in 'trap' 
units: m = u z = K = ks = 1. 



Eq. © is iteratively solved to self-consistency using the 
modified Broyden's method [3j. We work in a canoni- 
cal formalism which keeps TV and the total polarization 
P = (Nf — N±)/N fixed through the number equations 

N<r = f dzpa-(z). 



TEBD calculations 

To implement TEBD formalism, we employ a ID 
Fermi-Hubbard Hamiltonian to approximate the contin- 
uum quasi-lD polarized Fermi gases in harmonic traps: 



H = — J Y^2(4,<r c i-i,° + h - c -) + Uj2 n ht n hi 

a i=2 i=l 

L 



(7) 



BdG Calculation 

We treat H within the mean-field Bogoliubov-de 
Gennes (BdG) approach for which there are many excel- 
lent references [2j. Here we simply state the BdG equa- 
tions for the pair wave functions Uj(z) and Vj(z) which 
decouple H : 



Mt A ( z ) 
A(z) -Hf - 







Uj 




Uj 


H _ 








V 3 



(5) 



where Ej is the associated energy. Despite this, the 
BdG treatment has been shown to yield qualitatively re- 
liable answers [2j . In accordance with Fermionic commu- 
tation relations, the quasi-particle amplitudes are nor- 
malized as: J dz \uj (z)\ 2 + \vj (z)\ 2 — 1. In terms of which 
the Gap A(z) and the free energy Q, may be written as : 



A(z) = UT i u j (z)v* j (z)f(E j ), 



2=1 



(G) 



where/ (E) represents the Fermi-Dirac distribution func- 
tion: f(E) = l/(e E/kBT + 1). We follow a convention 
that N t > N h we define = ^/2// n and the FFLO 
wave number by qo = k F — k F . 

Our theoretical framework is encapsulated within 
Eqs. (JSJ) and © and we discretize the system of Eq. §5§ 
using a piece-wise linear finite element basis which en- 
sures the continuity of both u(z) and v(z). A reduction of 
Eq. © into even and odd parity states about z = is pos- 
sible due to anticipated reflection symmetry of A about 
this axis. Nevertheless, each independent sub-block with 
distinct parity presents a very large eigenvalue problem 
because of the slow convergence of Eq. (J6j) . The slow 
convergence is tackled using a hybrid BdG-semi-classical 
strategy similar to Ref. [2j . Starting from an initial state 



where L is the number of discretized lattice sites, c] a , c i cr 
are respectively the creation and annihilation operators 
for spin a particles at ith lattice site, J is the hopping 
amplitude between the neighboring sites, and U is the on- 
site interaction strength between two unlike spins. The 
connection between the Fermi-Hubbard Hamiltonan (|7| 
and the Hamiltonian (j4)) upon which the BdG calcula- 
tion is based can be seen as follows: In the trap units 

T 2 

we mentioned above, the hopping amplitude J = jrp, 
where I is the total length of system (in our dimension- 
less units). From these, the parameters -j = 2gxol/L 
and ^j- = — -j) 2 are choosen accordingly. In our 

calculation, typically we choose L — 300 ~ 400. With 
these properly choosen characteristic parameters, the dis- 
cretized Hubbard Hamitonian can be trusted to represent 
a continuum system as have been previously shown 

The TEBD algorithm utilizes the Schmidt decompo- 
sition and the convergence of the simulation is mainly 
controlled by the so-called Schmidt rank \, which is the 
number of eigenvalues retained when truncating the Hib- 
ert space. In our calculation, since the computation time 
scales as x 3 , the optimal value of x ~ 100 — 150 is chosen 
to ensure the convergence is good enough when compar- 
ing with the results with higher x- Another souce of error 
comes from the Trotter-Suzuki expansion for the time 
evolution operator. To reduce it, we adopt fifth order 
Trotter-Suzuki expansion in our calculation while choos- 
ing a small enough time step based on self-consistent sta- 
bility test. 
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